%% Michael D. Konig, Andrei Levchenko, Tim Rogers and Fabrizio Zilibotti
%% Aggregate fluctuations in adaptive production networks
%% Forthcoming in Proceedings of the National Academy of Sciences, 2022
 
%% Replication file (Matlab) for Figure 3, main text.

%% Requirements: 

% "sales_years_us_manuf.mat" -> Matlab data file with the following information (US manufacturing firms only): year, sales_firm_1,...,sales_firm_n
% "sales_years_us.mat" -> Matlab data file with the following information (US manufacturing firms only): year, sales_firm_1,...,sales_firm_n
% "sales_years_manuf.mat" -> Matlab data file with the following information (US manufacturing firms only): year, sales_firm_1,...,sales_firm_n
% "sales_years.mat" -> Matlab data file with the following information (US manufacturing firms only): year, sales_firm_1,...,sales_firm_n
% "in_deg_years_us_manuf.mat" -> Matlab data file with firms' in-degrees/number of suppliers (US manufacturing firms only).
% "in_deg_years_us.mat" -> Matlab data file with firms' in-degrees/number of suppliers (US firms only).
% "in_deg_years_manuf.mat" -> Matlab data file with firms' in-degrees/number of suppliers (manufacturing firms only).
% "in_deg_years.mat" -> Matlab data file with firms' in-degrees/number of suppliers (full sample).
% "out_deg_years_us_manuf.mat" -> Matlab data file with firms' out-degrees/number of customers (US manufacturing firms only).
% "out_deg_years_us.mat -> Matlab data file with firms' out-degrees/number of customers (US firms only).
% "out_deg_years_manuf.mat -> Matlab data file with firms' out-degrees/number of customers (manufacturing firms only).
% "out_deg_years.mat -> Matlab data file with firms' out-degrees/number of customers (full samlpe).
% "bonacich_years_us_manuf.mat" -> Matlab data file with firms' Katz-Bonacich centralities (US manufacturing firms only).
% "bonacich_years_us.mat -> Matlab data file with firms' Katz-Bonacich centralities (US firms only).
% "bonacich_years_manuf.mat -> Matlab data file with firms' Katz-Bonacich centralities (manufacturing firms only).
% "bonacich_years.mat -> Matlab data file with firms' Katz-Bonacich centralities (full samlpe).
% "lifetime_bonacich_log_likelihood.m" -> Log-likelihood function for lifetime.

clear all
close all

%% Parameters.
drop_non_us = true;
drop_non_manuf = true;
theta_0(1) = 0.1; % Initial value of delta (constant) for optimization algorithm.
theta_0(2) = 0.1; % Initial value of zeta (inverse shifted out-degree) for optimization algorithm.
theta_0(3) = 0.1; % Initial value of rho (indicator in-degree zero) for optimization algorithm.
lb = zeros(length(theta_0),1); %[]; % Lower bound for bounded optimization algorithm.
ub = [];  % Upper bound for bounded optimization algorithm.

opt_alg = 'fmincon'; % choose from: 'nm' 'fmincon' 'fminunc' 'sa'

%% Load sales data.
disp(['Load sale data from file...'])
if(drop_non_us)
    if(drop_non_manuf)
        load('./Data/sales_years_us_manuf.mat')
    else
        load('./Data/sales_years_us.mat')
    end
else
    if(drop_non_manuf)
        load('./Data/sales_years_manuf.mat')
    else
        load('./Data/sales_years.mat')
    end
end

sales_years_full = sales_years;
sales_years(:,1) = [];
nosales =  (sales_years==0) + (isnan(sales_years)==1);

N = size(sales_years,2);
startyear = zeros(N,1);
lastyear = zeros(N,1);
censored = zeros(N,1);
never_sales = zeros(N,1);
for i=1:size(nosales,2)
    if nosales(1,i) == 1
    %if no sales in 2003, start year is first year with sales
        if sum(nosales(:,i))<size(nosales,1) 
            startyear(i)=find(nosales(:,i)==0,1,'first')+2002;
        else
            %never any sales, just drop firm out of sample
            never_sales(i)=1;
        end
    else
         %otherwise, assume entry in 2003 (ignore left censoring)
        startyear(i) = 2003;
    end
    if nosales(end,i) == 0
        lastyear(i)=2002+size(sales_years,1);
        censored(i) = 1;
    else
        if never_sales(i)==0
            lastyear(i)=find(nosales(:,i)==0,1,'last')+2002;
        end
    end
end

%% Plot the empirical age distribution.
age = lastyear-startyear;
ind = find(age<=0);
age(ind) = NaN;
disp(['The number of observations for the empirical life time is ' num2str(length(age(~isnan(age))))])
disp(['Average empirical life time is ' num2str(nanmean(age))])
disp(['Standard deviation of the empirical life time is ' num2str(nanstd(age))])
disp(['Variance of the empirical life time is ' num2str(nanstd(age)^2)])
figure()
set(gca,'Layer','top');
set(gca,'FontSize',20);
set(0,'defaulttextinterpreter','latex')
box on
%h = histogram(y);
h = histfit(age,20,'exponential');
hold on
set(h(1),'facecolor',[0.5 0.5 0.5]); 
set(h(2),'color','k');
vline([nanmean(age) nanmean(age)],{'r','r'});
hold off
xlabel(['life time'])
title(['mean = ' num2str(nanmean(age))])
set(gca,'FontSize',20);
axis tight
if(drop_non_us)
    if(drop_non_manuf)
        print('-depsc','Figures/entry_exit_times_distr_us_manuf.eps')
        save('Figures/entry_exit_times_distr_us_manuf.fig')
    else
        print('-depsc','Figures/entry_exit_times_distr_us.eps')
        save('Figures/entry_exit_times_distr_us.fig')
    end
else
    if(drop_non_manuf)
        print('-depsc','Figures/entry_exit_times_distr_manuf.eps')
        save('Figures/entry_exit_times_distr_manuf.fig')
    else
        print('-depsc','Figures/entry_exit_times_distr.eps')
        save('Figures/entry_exit_times_distr.fig')
    end
end

%% Load covariates data.
if(drop_non_us)
    if(drop_non_manuf)
        load('./Data/in_deg_years_us_manuf.mat')
        load('./Data/out_deg_years_us_manuf.mat')
        load('./Data/bonacich_years_us_manuf.mat')
    else
        load('./Data/in_deg_years_us.mat')
        load('./Data/out_deg_years_us.mat')
        load('./Data/bonacich_years_us.mat')
    end
else
    if(drop_non_manuf)
        load('./Data/in_deg_years_manuf.mat')
        load('./Data/out_deg_years_manuf.mat')
        load('./Data/bonacich_years_manuf.mat')
    else
        load('./Data/in_deg_years.mat')
        load('./Data/out_deg_years.mat')
        load('./Data/bonacich_years.mat')
    end
end
in_deg_years(:,1) = []; % Remove the first column that contains the years.
out_deg_years(:,1) = [];
bonacich_years(:,1) = [];

%% Drop firms that have never any sales
todrop = (never_sales==1);

in_deg_years_matched = in_deg_years(:,never_sales==0);
out_deg_years_matched = out_deg_years(:,never_sales==0);
bonacich_years_matched = bonacich_years(:,never_sales==0);
sales_years = sales_years(:,never_sales==0);
startyear(todrop) = [];
lastyear(todrop) = [];
censored(todrop) = [];
n=length(censored);

%% Compute maximum likelihood life-time estimates.
disp(['Compute maximum likelihood life-time estimates...'])
switch opt_alg

    case 'fmincon'
        
        disp(['using a constrained optimization algorithm (fmincon)...'])
        [theta_opt,fval,exitflag,output,lambda,grad,hessian] = fmincon(@(theta) -lifetime_bonacich_log_likelihood(theta,n,in_deg_years_matched,out_deg_years_matched,bonacich_years_matched,startyear,lastyear,censored),theta_0,[],[],[],[],lb,ub,[]);
        
    case 'fminunc'
        
        disp(['using an unconstrained optimization algorithm (fminunc)...'])
        [theta_opt,fval,exitflag,output,grad,hessian] = fminunc(@(theta) -lifetime_bonacich_log_likelihood(theta,n,in_deg_years_matched,out_deg_years_matched,bonacich_years_matched,startyear,lastyear,censored),theta_0);
        
    case 'nm'
        
        disp(['using an unconstrained Nelder-Mead algorithm...'])
        [theta_opt,fval] = fminsearch(@(theta) -lifetime_bonacich_log_likelihood(theta,n,in_deg_years_matched,out_deg_years_matched,bonacich_years_matched,startyear,lastyear,censored),theta_0);
        
    case 'sa'
        
        disp(['using a constrained simulated annealing algorithm...'])
        [theta_opt,fval,exitFlag,out] = simulannealbnd(@(theta) -lifetime_bonacich_log_likelihood(theta,n,in_deg_years_matched,out_deg_years_matched,bonacich_years_matched,startyear,lastyear,censored),theta_0,lb,ub);
        
    otherwise
        
        error('No optimization algorithm specified.')
end

disp(['Parameter estimates for delta, zeta and rho: ' num2str(theta_opt)])

if(strcmp(opt_alg,'fmincon') || strcmp(opt_alg,'fminunc'))

    err = sqrt(diag(inv(hessian)));
    disp(['Estimated standard errors for delta, zeta and rho: ' num2str(err')])
    
    z = abs(theta_opt./err');
    pvalue = 2*(1 - normcdf(z));
    disp(['P-values for delta, zeta and rho: ' num2str(pvalue)])

    %% Plot rho estimates.
    if(length(theta_opt)>3)
        figure()
        set(gca,'Layer','top');
        set(gca,'FontSize',20);
        set(0,'defaulttextinterpreter','latex')
        box on
        rho_mean = theta_opt(3:end);
        rho_std = err(3:end);
        plot([0:1:length(rho_mean)-1],rho_mean,'-ok','LineWidth',1)
%         hold on
%         plot([0:1:length(rho_mean)-1],rho_mean+rho_std,'--r','LineWidth',1.5)
%         hold on
%         plot([0:1:length(rho_mean)-1],rho_mean-rho_std,'--r','LineWidth',1.5)
        hold off
        ylabel(['$\widehat{\rho}$'])
        xlabel('$d^-$')
        axis tight
        set(gca,'FontSize',20);
        if(drop_non_us)
            if(drop_non_manuf)
                print('-depsc','Figures/rho_estimates_us_manuf.eps')
                save('Figures/rho_estimates_us_manuf.fig')
            else
                print('-depsc','Figures/rho_estimates_us.eps')
                save('Figures/rho_estimates_us.fig')
            end
        else
            if(drop_non_manuf)
                print('-depsc','Figures/rho_estimates_manuf.eps')
                save('Figures/rho_estimates_manuf.fig')
            else
                print('-depsc','Figures/rho_estimates.eps')
                save('Figures/rho_estimates.fig')
            end
        end

    end
    
end
